In this tutorial we first solve the problem
\begin{cases} -u'' = f & \text{in }\Omega = (0, 1),\\ u = 0 & \text{on }\Gamma = \{0, 1\} \end{cases}
using standard (non block) dolfinx code.
Then we use block support of dolfinx to solve the system
\begin{cases} - w_1'' - 2 w_2'' = 3 f & \text{in }\Omega,\\ - 3 w_1'' - 4 w_2'' = 7 f & \text{in }\Omega \end{cases}
subject to
\begin{cases} w_1 = 0 & \text{on }\Gamma,\\ w_2 = 0 & \text{on }\Gamma \end{cases}
By construction the solution of the system is $$(w_1, w_2) = (u, u)$$
We then compare the obtained solutions. This tutorial serves as a reminder on how to solve a linear problem in dolfinx, and introduces multiphenicsx.fem.petsc.BlockVecSubVectorWrapper.
import typing
import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.mesh
import mpi4py.MPI
import numpy as np
import numpy.typing
import petsc4py.PETSc
import ufl
import viskex
import multiphenicsx.fem
import multiphenicsx.fem.petsc
mesh = dolfinx.mesh.create_unit_interval(mpi4py.MPI.COMM_WORLD, 32)
def left(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[np.bool_]:
"""Identify the left boundary of the domain."""
return abs(x[0] - 0.) < np.finfo(float).eps # type: ignore[no-any-return]
def right(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[np.bool_]:
"""Identify the right boundary of the domain."""
return abs(x[0] - 1.) < np.finfo(float).eps # type: ignore[no-any-return]
left_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim - 1, left)
right_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim - 1, right)
boundary_facets = np.hstack((left_facets, right_facets))
x0 = ufl.SpatialCoordinate(mesh)[0]
viskex.dolfinx.plot_mesh(mesh)
def run_standard() -> dolfinx.fem.Function:
"""Run the standard case: solve for the unknown u."""
# Define a function space
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 2))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
# Define problems forms
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
f = ufl.inner(100 * ufl.sin(20 * x0), v) * ufl.dx
# Define boundary conditions
zero = petsc4py.PETSc.ScalarType(0)
bdofs_V = dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim - 1, boundary_facets)
bc = dolfinx.fem.dirichletbc(zero, bdofs_V, V)
# Solve the linear system
u = dolfinx.fem.Function(V)
problem = dolfinx.fem.petsc.LinearProblem(
a, f, bcs=[bc], u=u,
petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})
problem.solve()
u.vector.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
# Return the solution
return u # type: ignore[no-any-return]
u = run_standard()
viskex.dolfinx.plot_scalar_field(u, "u")
def run_block() -> typing.Tuple[dolfinx.fem.Function, dolfinx.fem.Function]:
"""Run the block case: solve for the unknowns (w_1, w_2)."""
# Define a block function space
V1 = dolfinx.fem.functionspace(mesh, ("Lagrange", 2))
V2 = dolfinx.fem.functionspace(mesh, ("Lagrange", 2))
(u1, u2) = (ufl.TrialFunction(V1), ufl.TrialFunction(V2))
(v1, v2) = (ufl.TestFunction(V1), ufl.TestFunction(V2))
# Define problem block forms
aa = [[1 * ufl.inner(ufl.grad(u1), ufl.grad(v1)) * ufl.dx, 2 * ufl.inner(ufl.grad(u2), ufl.grad(v1)) * ufl.dx],
[3 * ufl.inner(ufl.grad(u1), ufl.grad(v2)) * ufl.dx, 4 * ufl.inner(ufl.grad(u2), ufl.grad(v2)) * ufl.dx]]
ff = [ufl.inner(300 * ufl.sin(20 * x0), v1) * ufl.dx,
ufl.inner(700 * ufl.sin(20 * x0), v2) * ufl.dx]
aa_cpp = dolfinx.fem.form(aa)
ff_cpp = dolfinx.fem.form(ff)
# Define block boundary conditions
zero = petsc4py.PETSc.ScalarType(0)
bdofs_V1 = dolfinx.fem.locate_dofs_topological(V1, mesh.topology.dim - 1, boundary_facets)
bdofs_V2 = dolfinx.fem.locate_dofs_topological(V2, mesh.topology.dim - 1, boundary_facets)
bc1 = dolfinx.fem.dirichletbc(zero, bdofs_V1, V1)
bc2 = dolfinx.fem.dirichletbc(zero, bdofs_V2, V2)
bcs = [bc1, bc2]
# Assemble the block linear system
AA = dolfinx.fem.petsc.assemble_matrix_block(aa_cpp, bcs)
AA.assemble()
FF = dolfinx.fem.petsc.assemble_vector_block(ff_cpp, aa_cpp, bcs)
# Solve the block linear system
uu = dolfinx.fem.petsc.create_vector_block(ff_cpp)
ksp = petsc4py.PETSc.KSP()
ksp.create(mesh.comm)
ksp.setOperators(AA)
ksp.setType("preonly")
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("mumps")
ksp.setFromOptions()
ksp.solve(FF, uu)
uu.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
ksp.destroy()
# Split the block solution in components
u1u2 = (dolfinx.fem.Function(V1), dolfinx.fem.Function(V2))
with multiphenicsx.fem.petsc.BlockVecSubVectorWrapper(uu, [V1.dofmap, V2.dofmap]) as uu_wrapper:
for u_wrapper_local, component in zip(uu_wrapper, u1u2):
with component.vector.localForm() as component_local:
component_local[:] = u_wrapper_local
# Return the block solution components
return u1u2
u1, u2 = run_block()
viskex.dolfinx.plot_scalar_field(u1, "u1")
viskex.dolfinx.plot_scalar_field(u2, "u2")
def compute_errors(
u1: dolfinx.fem.Function, u2: dolfinx.fem.Function, uu1: dolfinx.fem.Function, uu2: dolfinx.fem.Function
) -> None:
"""Compute errors between standard and standard block cases."""
u_1_norm = np.sqrt(mesh.comm.allreduce(
dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(ufl.grad(u1), ufl.grad(u1)) * ufl.dx)),
op=mpi4py.MPI.SUM))
u_2_norm = np.sqrt(mesh.comm.allreduce(
dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(ufl.grad(u2), ufl.grad(u2)) * ufl.dx)),
op=mpi4py.MPI.SUM))
err_1_norm = np.sqrt(mesh.comm.allreduce(
dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(ufl.grad(u1 - uu1), ufl.grad(u1 - uu1)) * ufl.dx)),
op=mpi4py.MPI.SUM))
err_2_norm = np.sqrt(mesh.comm.allreduce(
dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(ufl.grad(u2 - uu2), ufl.grad(u2 - uu2)) * ufl.dx)),
op=mpi4py.MPI.SUM))
print(" Relative error for first component is equal to", err_1_norm / u_1_norm)
print(" Relative error for second component is equal to", err_2_norm / u_2_norm)
assert np.isclose(err_1_norm / u_1_norm, 0., atol=1.e-10)
assert np.isclose(err_2_norm / u_2_norm, 0., atol=1.e-10)
print("Computing errors between standard and block")
compute_errors(u, u, u1, u2)
Computing errors between standard and block
Relative error for first component is equal to 4.923782490298636e-15 Relative error for second component is equal to 2.7591822783167812e-14